A numerical finite size scaling approach to many-body localization 
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We develop a numerical technique to study Anderson localization in interacting electronic systems. 
The ground state of the disordered system is calculated with quantum Monte-Carlo simulations while 
the localization properties are extracted from the "Thouless conductance" g, i.e. the curvature of 
the energy with respect to an Aharonov-Bohm flux. We apply our method to polarized electrons 
in a two dimensional system of size L. We recover the well known universal /3(g) = dlogg/dlogL 
one parameter scaling function without interaction. Upon switching on the interaction, we find that 
/3(g) is unchanged while the system flows toward the insulating limit. We conclude that polarized 
electrons in two dimensions stay in an insulating state in the presence of weak to moderate electron- 
electron correlations. 
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Since the early days of Anderson localization it is 
believed that in the thermodynamic limit, an arbitrary 
small disorder is enough to drive a two dimensional elec- 
tron gas toward an insulator [2]. At the origin of this 
prediction is the scaling theory of localization 0] which 
conjectured that the evolution of the conductance with 
the system size obeyed a simple one parameter scaling 
function. An important numerical effort has since been 
devoted to establish the presence of this scaling 0, H| and 
calculate the scaling function. While this one electron lo- 
calization picture is now reasonably well understood, the 
corresponding many-body problem, where not only dis- 
order but also electron-electron interactions are consid- 
ered, is yet unsolved. An important litterature has been 
devoted to the very strong [5j and weak disorder limit 
[si 0] but very little on the interplay between interac- 
tion and localization itself Q. It was generally assumed 
that electron-electron interactions did not modify dras- 
tically the one electron physics, so that the observation 
of a metallic state in two-dimensional Si MOFSETs in 
1994 [9(] came as an important surprise. It gave rise to 
a new interest in the subject 13, U, 12 1 and raised the 
question of the possibility that electron-electron interac- 
tion could stabilize a metallic phase. Recent progresses 
in the weak disorder limit seem to indicate that it could 
indeed be the case[13| . 

Numerical methods have proved to be very useful 
in putting the scaling theory of localization for non- 
interacting particles on very firm grounds. It is therefore 
very tempting to try to develop similar approaches for 
the many-body problem in spite of the intrinsic difficul- 
ties in dealing with correlations. Indeed, a number of 
technical problems need to be overcome, (i) Obtaining 
the ground state of a decently large number N of corre- 
lated particles is already a challenging task, (ii) In order 
to study Anderson localization (and not a mere trapping 
of the electrons which is found for very strong disorder) 
one needs the localization length £ to be rather large 
yet smaller than the system size L which must hence be 
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FIG. 1: (Color online). Scaling function f3(g) as a function 
of logg, in rectangular systems with v — 1/24. The sym- 
bols correspond to different system sizes, TV = 16 particles in 
16x24 sites (empty symbols), TV = 25 in 20x30 (full symbols), 
N = 36 in 24 x 36 (striped symbols) , and different strengths of 
the interaction, r s = (circles), r s — 2 (squares), r s = 4 (up 
triangles) , r s — 6 (diamonds) , for various strengths of the dis- 
order (0 < r w < 30). The black lines are the expected asymp- 
totic limits. Inset: idem for various filling factors, v — 1/24 
(circles), v = 1/54 (stars) and v — 1/96 (triangles down), 
with TV = 16 particles, at r s = 0. Upon increasing disorder 
or interaction, the system flows toward the insulating limit as 
indicated by the arrow. 



rather large itself, (iii) One needs to calculate a physical 
observable sensitive to localization which must hence be 
some sort of correlation function [3]. Indeed, thermo- 
dynamic quantities such as the electronic density do not 
show localization in average. 

In this letter, we propose a practical scheme to study 
numerically many-body localization. We use zero tem- 
perature Green Function quantum Monte-Carlo (GFMC) 
technique to study the ground state of the system. This 
technique takes full advantage of the fact that the non- 
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FIG. 2: (Color online), logg as a function of the system 
size L, in square systems with ^=1/25. Disorder and inter- 
actions are increased from top to bottom as indicated by the 
arrow. The different curves are: r m =0.005, r s =0 (x), r w =10, 
r s =0 (+), r w =15, r s =0 (*), r w =l7.5, r s =0 (0), r w =20, r s =0 
(A), r w =2Q, r a =2 (a), r w =20, r s =4 (►), r w =20, r s =6 (T), 
r„=22.5, r s =0 (□), r m =22.5, r s =2 (■), r™=22.5, r s =4 (♦), 
r w =27.5, r s —0 (o) and r TO =27.5, r 3 —2 (•). The collapse of 
the linear fits (dashed lines) at L = is the signature of one 
parameter scaling in the localized regime. 



interacting ground state can be found exactly (by di- 
agonalization of the one-body problem) so that upon 
switching on the interaction the GFMC simulations are 
done with a very good starting point. Our tool to mea- 
sure the localization properties is the "Thouless conduc- 
tance" of the system, which can be related to the dis- 
tribution of the winding numbers in the imaginary time 
path integral. We apply our method to polarized (spin- 
less) electrons in two dimensions for which both the- 
ory and experiments agree that the system is insulating. 
The main point of scaling theory of localization is that 
P(g) = dlogg/dlogL which depends on disorder, inter- 
action, density and size is in fact a function of g only. 
Our chief result is presented in Fig. [T] where we establish 
this scaling for interacting electrons. We find that (3(g) 
is unaffected by the presence of the correlations due to 
Coulomb repulsion in agreement with what is expected 
from the weak disorder limit [IE] ]. Upon increasing the 
interaction strength, the system flows toward the insu- 
lating limit and the system localization length £ (shown 
in Figd]) decreases. 

Model and method. We consider a system of TV spinless 
electrons in a rectangular L x x L v lattice with periodic 
boundary conditions. The Hamiltonian reads, 

U 



restricted to nearest neighbors and the den- 

sity operator. The disorder potential vp is uniformly dis- 
tributed inside [-W/2, W/2]. U is the effective strength 
of the two body interaction V(r). To reduce finite size ef- 
fects, V(r) is obtained from the bare Coulomb interaction 
using the Ewald summation technique. The expressions 
for V(r) has been given in [16j |. At small filling factor 
v = N/(L x L y ) <C 1, we recover the continuum limit 
and we are left with two dimensionless parameters, the 
usual r s = m* e 2 / (h 2 ey/wri) (m* effective mass, e elec- 
tron charge, e dielectric constant and n electronic den- 
sity) interaction parameter which for our model reads 
r s = U I (2ty/TTv) and a parameter r w = W/ (t^/v) con- 
trolling the strength of the disorder. In the diffusive limit 
without interaction, the product of Fermi momentum kp 
by the mean free path I is given by kpl = 1927r/r^. 

The GFMC method and our particular implementa- 
tion has been given in [161 ] to which we refer for de- 
tails and references. GFMC is a lattice version of the 
standard zero-temperature quantum Monte-Carlo meth- 
ods (like diffusive quantum Monte-Carlo) that have en- 
joyed important success for both bosonic and fermionic 
systems [171 ]. Its principle is to project an initial varia- 
tional guiding wave- function (GWF) \^g) onto the ex- 
act ground state \^o) by applying the projector operator 
e -pH m a stochastic way. Quantum Monte-Carlo meth- 
ods suffer from the so called sign problem when dealing 
with fermionic statistics. One way out of the sign prob- 
lem which has been quite successful is the fixed node 
approximation 17] where upon projection onto I'J'o) the 
sign of the wave-function is kept fixed. The method is 
variational and calculates the best wave-function com- 
patible with the nodal structure of the GWF. Impor- 
tant effort is usually spent looking for a GWF as close 
to the real ground state as possible. In the present case 
however, we can obtain the ground state without interac- 
tion exactly by diagonalizing the corresponding one-body 
problem. As we are interested in the evolution of the lo- 
calization properties upon switching on the interaction, 
we have an excellent starting point to begin with. The 
general form of our GWF is a Slater determinant multi- 
plied by a Jastrow function, 

y G {r u r 2 ...r N ) = Det X J[ J(\n - fj\). (2) 

i<j 

The Jastrow part introduces some correlation and ac- 
count for Coulomb repulsion. We use modified Yukawa 



J(r) 



exp[ 



aA(r s 



(1 



-B(r B )r/a 



)] 



functions [18 

where a = 1/ ^Jtu) is the average distance between elec- 
trons. A(r s ) and B(r s ) are variational parameters that 



H = —t ^2 c \ c r' + vpnr + — ^ V(r — r')npnp + A, we optimize while imposing the cusp condition B 



(r,r') 



(1) 

where ct et cp are the usual creation and annihilation 



\Jr s /A to reproduce the short distance behaviour. The 
nodal structure of the GWF depends only on the Slater 
determinant Det [<^i(r*j)] which enforces the antisymme- 
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try. The Slater determinant is constructed out of one- 
body orbitals 4>i that are obtained in two different ways 
leading respectively to \I/ii q and ^Har GWF. The orbitals 
of \l/ii q are calculated by exact diagonalization of the one- 
body (disordered) problem so that \I/iiq coincides whith 
the exact ground state without interaction (r s = 0) . The 
calculation of the orbitals of ^Har proceeds in a similar 
way but wc include iteratively the (Hartrcc) mean field 
potential due to the density (n?) of electrons in the one- 
body problem. The Hartree potential tends to screen the 
disorder leading to an increase of the GWF's localization 
length. We shall verify however that the GFMC results 
are not sensitive to the choice of GWF. 

Measuring the localization properties. The idea to use 
the sensibility D of the system to a tilt in its bound- 
ary conditions as a criteria of localization was introduced 
very early by Edwards and Thouless [19(. Indeed, in a pe- 
riodic system, the position of the boundary can be moved 
by a simple gauge transformation so that all sites are 
equivalent with respect to the boundary and a localized 
state is expected to have an exponentially small sensitiv- 
ity to the boundary. More precisely, in presence of a small 
Aharonov-Bohm flux <f>, a current I = —dE/d<f> flows in 
the system (E is the total energy). When the flux is 
small, we have I = —Dcf> where D = d 2 E/d<fi 2 \ < p = o is the 
curvature of the energy. This quantity D is referred as 
the "Thouless conductance", the Drude weight, the con- 
ductivity stiffness or the superfluid stiffness depending on 
the context. For bosonic systems, D is simply related to 
the superfluid fraction [13]. For fermions, it is related to 
the low frequency limit of the imaginary part of the con- 
ductivity [2l(. For disordered system the product of D 
with the density of states is proportional to the conduc- 
tance of the system 13, 22 1. In most cases, D is positive 
and the system is diamagnetic but in some instances (one 
dimensional systems with even number of spinless elec- 
trons [23l | . or two dimensional systems with degenerate 
ground states (24j) a paramagnetic (D < 0) response can 
been found so that the widely used interpretation of D 
as a conductance can sometimes be problematic. In any- 
case, it is a good measure of the localization properties 
of the system. 

Following [23], we calculate the diffusive constant 
g of the motion of the center of mass of the sys- 
tem in imaginary time along the x direction, g = 
lim^i N{Rl((3))/(t(3), where (R 2 x {(i)) is the second mo- 
ment of the center of mass along x. An example of the 
calculation of g is shown in the left panel of Fig|3] g is 
simply related to D as [13] g = DL 2 /(Nt) yet g is easier 
to access in the simulations. We note that in the fixed 
node approximation, g is always positive by construction. 
In particular, in the absence of disorder r w = and inter- 
action r s — 0, we find 3 = 2, i.e. the sum of the curvature 
of the individual one-body levels 21 1. 

Numerical results: square samples. We now turn to 
the numerics and show that g is an appropriate mea- 
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FIG. 3: (Color online). Left: example of (R%.(/3)} for four 
different samples, r w = 20 and r s = 4. The thick lines are 
the fits used to extract the values of g. Upper right: f/a 
as a function of r w , for v — 1/25 (circles) and v = 1/54 
(squares) at r s = 0. Lower right: £ calculated with GFMC as 
a function of the exact result £ ex at r s = for ^ = 1/25 and 
20 < r w < 30. The dashed line is a linear fit. 



sure of localization. In Fig. [2~] we plot log g (averaged 
upon disorder) as a function of system size L for square 
samples at fixed filling factor. Without interaction, scal- 
ing theory of localization [2]] predicts that log g is inde- 
pendent of L for small disorder (Ohm's law) while for 
strong disorder, g decreases exponentially with L so that 
log g = log go — L/t;. The existence of a universal scaling 
law 0(g) in this case means that log go is just a constant, 
independent of the disorder strength. Indeed, the nu- 
merics are fully consistent with this picture: for r w < 15, 
log g is roughly constant. Upon increasing disorder fur- 
ther, log g starts to decrease linearly with L while all 
curves intercept at a single point at L = 0. We find 
loggo = — 0.47±0.08. Further check of the method can be 
done by comparing the localization length £ with results 
^ ex of an exact diagonalization of the one-body Hamil- 
tonian [25|. The result is shown in the lower right panel 
of Fig. 03 We find both methods in good agreement as 
£ oc (£ ex — 1) [13] • In the upper right panel of Fig. [3l we 
plot £/a for various values of W and v and verify that it 
is a function of r w — W/(ty/v). 

We are now ready to switch on the interaction in our 
system. We find (Fig. [2~1 full symbols) that upon in- 
creasing r s , the localization length decreases so that the 
system becomes more insulating. More importantly, all 
the curves still intercept at the same single point at L = 
indicating that the universal scaling function (3(g) is un- 
affected by the electron-electron correlations. To confirm 
this important point, we have performed simulations with 
our two different wave-functions ^fn q and ^Har and find 
that the resulting £ agree very well as shown in the low- 
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FIG. 4: (Color online). £/a as a function of r s for various 
disorders, r w = 19.6 (circles), r m = 22 (diamonds) and r w = 
24.5 (triangles). Top: N = 25 in 20 x 30 sites. Bottom: 
N — 16 in 16 x 24. Empty (full) symbols correspond to ^Har 

(*Uq). 



find that the universal scaling function [3(g) is unaffected 
by the interactions. Yet, upon increasing the interaction 
strength, the system flows toward the insulating limit. 
This picture is in agreement with what is expected in the 
weak disorder and weak interaction limit [la ] . A natural 
extension of this work would be the study of non polar- 
ized electrons where the existence of an intrinsic metal- 
insulator transition remains a controversial issue. We 
note that in Fig. [4j the localization length £ extrapo- 
lates to zero at r s as 10 for the two studied values of 
disorder. This could be the signature of a transition to- 
ward some sort of disordered Wigner crystal. Remark- 
ably, r s w 10 also corresponds to the density at which 
the metal-insulator transition was observed for all but 
the cleanest samples [27j . 
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est panel of Fig. |U This is a good test of the robustness 
of the method: even with ^Har who tends to be less lo- 
calized than Wiiq, the FN-GFMC algorithm succeeds to 
find the (more localized) ground state. 

Numerical results: rectangular samples. The scaling 
function (3(g) = dlogg/dlog-L can in principle be ex- 
tracted from Fig. [5] by finite differences. More pre- 
cise results are obtained using rectangular (L y = l.5L x ) 
samples. We now calculate the diffusion constants g x 
and g y along the two different directions and compute 
(3(g) = \o g(g y /g x )/\og(L v /L x ) as a function of logg = 
log(g x g y ) /2. This scheme allows us to obtain the full 
scaling curve for a single system by varying the disorder 
parameter r w . The result is shown in Fig. [1] for vari- 
ous values of N, r w and r s while different values of v are 
shown in the inset. All data collapse on one single curve. 
We emphasize that no operation is needed to obtain this 
collapse, Fig. [T] shows raw data. The asymptotic curves 
are simple straight lines of slope one and zero which in- 
tercept at log go (which has been extracted from Fig. [2}. 
Small deviations from scaling is observed for the small- 
est size N = 16 at r s > 6. Much larger deviations were 
found for N — 9 particles (not shown). For N = 25 
and higher, the collapse was perfect up to our statisti- 
cal accuracy. The corresponding localization lengths are 
plotted in Fig. |4j They decrease nearly linearly with r s , 
up to r s — 7. Although it is very difficult to reach higher 
values of r s , it is very likely that the localization length 
stays below the non interacting one which rules out the 
possibility of a metallic behaviour. 

Conclusion. We have proposed a practical method to 
study Anderson localization in presence of many-body 
correlations. For polarized two-dimensional electrons we 
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